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A new concept called here the "trace function" was formu- 
lated and developed in this work. 

The notion is that the phase differences/time delays of 
the signals received from an array of sensors which has a given 
defined geometry describe a certain "trace function" which is 
constant for that array at a given look direction. If one can 
find a geometry for which the trace function's basic shape is 
not dependent on the look direction, but only some of its 
parameters are, then one has a powerful method to correlate 
(compare) the trace function of the received signal to a stored 
replica. 

The concept of the trace function was formalized and its 
application to some typical array geometries was demonstrated. 
Furthermore, it was shown that for highly symmetric arrays like 
the circular and linear arrays the trace function reduces to a 
particularly simple form. 

This characteristic was used to derive a relatively simple 
and manageable MMSE estimator for the bearing of the incoming 
signal. The estimator is applicable to either narrow band 
Signals by use of phase difference trace functions or toa 
wide band signal using the time delay trace function. 

The performance of the estimator was checked by simulation 
and compared to the CRLB as adapted to this application. 

Finally, a system configuration applying trace function 
principles was outlined and the major problem areas caused by 
the specific application were identified, reviewed and some 
suggestions to solutions were made. The principle can be 
applied to any other situation in which phase differences 


and/or time delays within an array of sensorscan be measured. 
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۳ PURPOSE 

The estimation of the bearing, frequency and amplitude 
of acoustic signals for passive sonar is normally accomplished 
With an array of sensors which has a defined geometrical 
arrangement. In most cases there are several factors which 
influence the geometry, of which the most prominent are the 
size, number of elements and frequency coverage, versus the 
desired performance. As a result of various tradeoffs, tra- 
ditionally three main geometrical configurations have been 
applied: 

~ the linear array 

~ the circular or spherical array 

- the conformal array. 
The linear array is typically applied in fixed sites and in 
towed configurations and as such the size (length) is not the 
dominant limitation. The circular or spherical array and the 
conformal array are applied normally in moving platforms and 
the size limitations (maximum aperture) is a major factor in 
the design. The size limitation implies a limitation of an- 
gular accuracy and resolution at low frequencies. The extent 
of the performance limitation can be reduced by proper pro- 
cessing of the array outputs. 

To achieve the desired performance, the processing of 


array data is done in two domains: space domain and time 





domain and their derived (by transform) versions; spatial 
and temporal frequency. 

Several approaches have been taken to the above pro- 
cessing, from the simplest one which is the delay and sum 
beamformer [35, 37, 33] through the correlation processor 
۳۳ 2۱۱ to the optimal processing l 1, 2 , 5 5 135 20> 245 39] 
with adaptive version as the "ultimate" processor [12, 20, u2]. 
All those techniques solve to some extent the processing pro- 
blem of arrays whose size is such that the interelement 
Spacing is on the order of A/2 for the lowest frequencies of 
interest, and have enough elements to sample adequately the 
space/time field. 

Attempts to overcome the spatial resolution problem have 
been made by applying the "superdirective array" concept. 
Although this technique gives infinitely good spatial resolu- 
tion for noise free environment, it has not been found to be 
practical because of poor aperture efficiency and sensitivity 
of the solution to tolerances in implementation [11, 35, 37]. 
Recently some constrained sensitivity solutions to the super- 
directive arrays were presented which give a less sensitive 
and more efficient array [10, 24]. 

Some of the above processors are implemented as "beam- 
formers" while the others are considered as estimators of the 
incoming signal direction. As a general conclusion it can be 
said that an adequate solution to the problem of small arrays 
at low frequencies has not yet been fully derived. The aim 


of this dissertation was to develop such a solution. 





B. MOTIVATION 

The practical situation which motivated the initiation 
of this work was the small circular and cylindrical arrays 
found in a great number of submarines and surface ships. 
These arrays are constrained in size to few meters diameter 
which in turn limits their ability to perform spatial pro- 
cessing of low frequency signals. The size of the whole 
array is small compared to the wavelength at frequencies of 
BOUSHZ or less. 

However, those low frequencies are the very ones which 
offer longer detection ranges and better classification means. 
The spatial processing problem can be solved by adding a long 
line array (towed in general) which will give the needed 
aperture for the low frequency spatial resolution. For small 
platforms this approach may not be practical and a different 
way has to be found by signal processing of the small circu- 


lar array. 


EOD BINE OF THE DISSERTATION 

A novel idea which is believed to be original is used as 
the basis for this dissertation. This idea is basically a 
new way to look at the available information from an array 
of sensors called the "trace function" by the author. 

It is based on the fact that the directional information 
of an incoming plane wave is mainly contained in the received 
Signal phases, or to be more precise, in the relative phase 
differences between received signals at the various elements 


of the array. 
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Those phase differences are representative of the time 
delay of the plane wave between the elements of the array. 
The notion is that the phase differences/time delays of the 
signals received from an array of sensors which has a given 
defined geometry describes a certain "trace function" which 
is constant for that array at a given look direction. If 
one can find a geometry for which the trace function's basic 
shape is not dependent on the look direction, but only some 
of its parameters are, then one has a powerful method to 
correlate (compare) the trace function of the received signal 
to a stored replica. 

Fortunately one such geometry is provided by the circular 
array, which was the motivation behind this work and whose 
phase difference trace function is a simple two dimensional 
Sine-cosine "wave" of exactly one period of each dimension. 

The apparent "phase" of this trace function is indicative 
of the bearing of the incoming plane wave signal. The trace 
function idea can be used with array geometries other than 
the circular one, however, it will be less efficient as the 
symmetry is broken. For example, it will be shown in later 
Chapters how it applies to linear arrays. 

General Definition - A trace function is a discrete plot 
of phase difference/time delay versus element indices for all 
element pairs. 

Fig. I-1 shows such a trace function for a random array 
while Figs. I-2 and (-3 show one each for a linear and 


circular array. 
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Two main features can be readily observed for any trace 
fone tion: 
- the elements on the diagonal are always zero (AT,,=0 
and Ap. .=0). 
- the trace function is always antisymmetric across the 


diagonal ( ar سر د‎ i 


and 7 2-4 ュ フ š 

For the noise free conditions the trace function can be 
derived easily for a given geometry. Trace functions are 
defined separately for each look direction and for each fre- 
quency resolution cell (or for wide band time delay). 

In Chapter II of this dissertation the subject of defini- 
tion and construction of the trace function will be treated 
in great detail. 

For the received signal in a noisy environment the above 
trace functions must be estimated for each narrow band signal 
present and for the wide band data. 

Once the trace function of the received signals is 
estimated it can be "matched" to a stored replica of the 
noise-free trace function of the array (either in discrete 
or analytic form). 

The task of estimating the trace function from the re- 
ceived data is not trivial mainly because it has to be done 
in the very difficult conditions of small array size and 
highly correlated noise between sensors. 

All optimal processors use some scheme to "whiten" the 
noise in spatial domain and/or temporal domain before it is 


further processed. This approach has to be considered in 
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this case too. However, the trace function concept poses a 
heavy demand on that type of preprocessing because the 
preservation of relative phase within the array is required. 
Present prewhitening schemes are not readily applicable and 
ways have to be found to perform the above task. Suggestions 
about relative phase preserving spatial decorrelation methods 
are made in Chapter V (System Considerations) which use a 
compensation scheme to preserve phase, and are based on the 
knowledge (a priori or estimated) of the ambient noise spec- 
trum and covariance functions. 

The estimation of phase differences can be done by es- 
tablished cross spectral estimation methods, especially the 
"short modified periodogram" method [40, 25, 26] which was 
adopted in this work and is reviewed in Chapter V for this 
application. 

For the time delay estimation several approaches are 
possible and this subject is treated extensively in the 
current literature [6, 9, 14, 15]. 

Of the above methods, the one which was the most appli- 
cable 1s reviewed in Chapter V and some adaptations to this 
work are suggested. 

The match of the estimated and stored trace function was 
done by a minimum mean square error (MMSE) estimator which was 
fully derived in this work. The derivation of the estimation 
is given in Chapter III and its performance is detailed in 
Chapter IV. The performance of the estimator was checked by 


computer simulations and the results are presented and compared 
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to a classical benchmark, the Cramer-Rao Lower Bound (CRLB). 

A block diagram for a possible implementation scheme is 
given in Chapter V. 

۱ VI presents a concise summary of this work. In 
addition it is suggested that the trace function concept and 
the estimators derived in this work are not limited in their 
application to sonar. They can be applied to related areas 
like seismic, and ocean wave analysis. The chapter concludes 
with some proposed topics for future work which were opened 


up by the introduction of the trace function principle. 


D. RELATED WORK 

As the concept is new, very little related work can be 
mentioned as applying strictly to the basic idea of the trace 
function. Only two works [3, 22] can be defined as such. 

The rest of the references used in this dissertation could be 
more correctly described as supporting works. 

Two previous papers published by the author (with G. L. 
Sackman) [30, 31] on the subject of the dissertation are also 
cited. 

The two paper's which are related to the subject [3, 22, 
op.cit.] were brought to the attention of the author recently 
after the trace function concept had already been developed. 
E c'at1ionshrp is restricted to the general concept of 
"matching" some function of the cross spectral matrix to an 


incoming plane wave. Those two papers were analyzed in order 
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to contrast with and emphasize the present work. 

Mee wtirst paper by Munk et al [22] is the basic one in 
which the authors address the problem of estimation of the 
direction of ocean waves arising from distant storms using a 
relatively small array of three sensors. The problem is 
somewhat similar to the present one, however it is ina 
totally different discipline and the methods used for solu- 
tion are also significantly different. The solution which 
was shown in that paper is a Meme - fitting of a full; 
measured cross spectral matrix (complex in general) of an 
array, to anassumed form of such a matrix for a plane wave. 
Their method used an "educated" search technique to minimize 
the error in the fitting process. 

Very good results were achieved, however the computa- 
tional complexity was tremendous. 

The second paper by Bennett[ 3] is a tutorial report 
published almost ten years after the first one. In this 
report the author formalizes and reviews the methods of the 
paper by Munk et al adding some alternatives and demonstrating 
the theory by the results of processing real data. However, 
no departure was made from the basic approach of Munk et al. 

It is worthwhile to mention that Bennett concluded that 
this method of least square fitting was the best available 
at that time for bearing estimation. In contrast to the 
above the present work uses only certain relevant parts of 
the cross spectral matrix - the phase differences - as a 


measure to be fitted and introduces the notion of the trace 
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function as a replica which is only dependent on array geo- 
metry. This simplification leads to an explicit MMSE 
estimator equation for some well defined array geometries. 

The option of "educated" search remained open, but for a 
much simpler set of data. All this is done without degrading 
performance; in fact it is demonstrated by simulation that 
the present estimator is asymptotically efficient (achieves 
C. R. Lower Bound). 

Another fact which should be mentioned is that the deri- 
vation of the trace function concept was initiated from the 
end product, the bearing, and how it effects the various 
measurable parameters in the array and not from the mathemati- 
cal structure of the cross spectral matrix as in the previous 
works. 

Nevertheless, the merit of these related works remains 
valid and both the previous and the author's methods comple- 
ment each other, opening perhaps a whole new way of pro- 
cessing array data. The supporting works will not be con- 
Sidered in this section mainly because they cover a wide 
variety of subjects and do not affect the heart of the pro- 
blem. They are mentioned and referenced elsewhere in this 


work, each along with its applicable subject. 


۳ CONTRIBUTIONS 

The principal contributions of this dissertation can be 
summarized as follows: 

1. The first introduction and derivations of the phase- 


difference/time-delay trace function concept for arrays of 
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TOES, sts applicability and limitations. 

2. Derivation of a MMSE bearing estimator which 
implements the trace function concept for circular and 
linear arrays, and its performance. 

3. Definition of a system which applies the new concepts 
along with suggestions for solution of some specific problems 
which arose from the special demands of phase preserving 


preprocessing for the trace function concept. 
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iad O ERENCE LIME DELAY TRACE 
۳ ۱۱۳ OF ARRAYS - THE CONCEPT 


A. INTRODUCTION 

The basic concept of the trace function of an array is 
the basis for the bearing estimation method developed in this 
dissertation. 

In order to give the reader a good understanding of the 
subject, the concept of the trace function is defined and 
developed for a general three dimensional array first and 
then for some typical two dimensional geometries of arrays. 
Finally the advantages of this concept will be pointed out 
for symmetrical arrays, especially the circular array. 

Shown in Fig. II-l is an arbitrary three-dimensional 
array of N receiving elements. The location of the i-th ele- 
ment with respect to a fixed (x,y,z) cartesian coordinate 
system is defined by the vector P. specifying both distance 
and direction from the center of the coordinate system. 

A uniform plane wave propagates with speed c in a 
direction a with respect to the coordinate system. The 
vector a has a unit length. 

The plane wave is assumed to contain a wide band random 
signal, Ss and a set of M line (sinusoidal) components with 
random initial amplitude and phase, S b: Then, 1f we assume 
that the same plane wave arrives at each element, we can 


Specify vector signal to be: 
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Uniform Plane Wave Target Signal 
Impinging On An Arbitrary 
Three Dimensional Array Of Elements 
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The above definition of a uniform plane wave target signal 
shows that the signal component at each element is identical 
vT ior a delay. This is under the assumption that atten- 
uation is negligible and no change in the speed of sound occurs 
within the array dimensions 

From a total system point of view we have to include also 
the additive noise component at each element; however, for 
the conceptual demonstration of the trace function the noise 
free condition is adequate. Under the above stated conditions 
it is suggested that these delays لر ج)‎ are the parameters 
which contain the directional information of the plane wave. 

It is convenient to consider the wide band portion of the 
signal and the line components separately, although it will 
be shown that the trace function concept is applicable to both 


in a Similar manner. 


B. NARROW BAND COMPONENT 
As stated before, the narrow band component S b 1s com- 


posed of a set of sinusoidal lines: 


S plti) = > A, COS Lw, (t-1,) + ل و‎ (II-3) 


M 
where Ay - line amplitude 
$5 - line phase 
The phase difference trace function is constructed by 
taking all possible pairs of elements and observing that the 


phase difference between any pair of elements' signals S prio 9) 


and Sp 59 is: 
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C = wy 0 (II-5) 


Using (II-2) and (II-4): 


Ad. (II-5) 


EE = TT.) = wy = 
where "i" and "j" are equivalent indices of the array ele- 
ments and denote the pairs which are considered for phase 
difference or time delays. 

Definition: A phase difference trace function is defined 
as the plot of Ad (phase difference as a discrete function of 
"i" and "j" (element numbers in the array) for "2" (frequency) 
constant. 

An example of such a trace function is shown in Figure 
II-2 for a random array. 

Observations: 

a. Each Ab is a constant for a given look direction 
a so that for a fixed geometry (fixed D.) and frequency w, 


the trace function in invariant in time. 


b. For various frequencies w,=1,2..M we get a family 


k 
of such invariant trace functions which are linearly related. 
c. For various look directions we get a set of such 
families, one for each look direction a. 
From the above observations we can conclude that if we 
know such a set of families of trace functions which are 


invariant for a given array geometry we can "match" a re- 


ceived and estimated trace function to the precomputed members 
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of the sets/families and find the look direction for each 
frequency resolution element. 

The EES of performing the above match can be quite for- 
midable for a practical situation, unless some regularity and 
symmetry can be found within the sets. It may even be pro- 
hibitive to be processed in real time systems. 

Fortunately, for some practical and symmetrical geometries 
of arrays, this task of matching is facilitated. 

It is easy to see that this concept applies to multiple 
targets as well if we assume that we can always find a fine 
enough frequency resolution such that two different targets 
will have their line components in different frequency anal- 
ysis resolution cells. (This condition is easily achieved 
in practice.) If they are coming from different look direc- 
tions they will be members of different families as defined 


above and they will be recognized as such. 


C. WIDE BAND COMPONENT 

For a wide band random process the phase is not defined, 
hence we will use the time delay as a measure for the trace 
function. As mentioned before the wide band component of the 
Signal is defined as Sp (tT). 


In analogy with the line components 


6 ° um) 
= e - لي لاگ رل ېي‎ a 
تت۸‎ (1: E > (II-6) 


1J 


Definition: A time delay trace function is defined as the 
plot of AT. (time delay) as a discrete function of "i" and 


"j" (element numbers in the array). 
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Figure II-3 shows such a time delay trace function for 
a random array. 

The observations for the narrow band components are 
readily applicable to the wide band components as well with 
one exception: the frequency domain collapses into only one 
trace function, such that it results ina single set of trace 
functions, one for each look direction a. 

The same conclusions as for the narrow band components 
apply, however the multiple target distinction based on 
frequency resolution cannot be done any more. Some sugges- 


tions to overcome this problem will be made in Chapter V. 


D. TRACE FUNCTION CALCULATION 
Mmenalytic Trace Functions 
The trace functions of arrays which have a geometry 
such that they can be expressed in a closed form equation are 
defined here as "analytic". The calculation of these trace 
functions is simple and is based on geometry only. In this 
work only two dimensional arrays were treated, however the 
extension to three dimensions is straightforward. 
The calculations are derived from equations (II-5) 
and (II-6) for phase difference and time delay trace functions 
respectively. 
These derivations are shown in Appendix A for circu- 
lar and linear arrays. 
a. Circular Array 
The configuration of this array is shown in Figure 


II-4. From Appendix A: 
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- Time delay trace function 


R š T T 5 T ووو‎ 
e e = سے ہے‎ — + = ーー = emm 
ATs. 22 sin ۳ (1+3) ١| هذه‎ [ (il D| (II-7) 


- Phase difference trace function 


_ ۹ Ge Dee 
Absa ig = ۵٢د٥‎ = PM ات‎ d 
WER Ee 
x sin [$ (43) | (II-8) 
where R,D - radius and diameter respectively 
N - number of elements in the array 
0 - incoming signal bearing (look direction) 





Figure II-4 


Circular Array Configuration 
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In Figure II-5 through II-13 are shown the three 
dimensional plots of the phase difference trace functions - 
for a circular array of 32 elements (Equation II-8) 


Ag 


at various look directions (bearings). The effect of fre- 


isis S 


quency change can be seen by comparing Figures II-5 and II-14. 
This effect is onlyonthe "amplitude" of the trace function 
but not on its basic form (similar effects are caused by 
changes in diameter-D or sound velocity-C). The computer 
program for ploting these trace functions on a Hewlett-Packard 
9845T computer is presented in Appendix B-1l. 

Observations: 

(1) The trace functions of circular arrays 
are always a sampled two dimensional combined cosine and sine 
surface of exactly one period (360°) in each dimension ("i" 
and "j"). 

(2) The "phase" of the above two dimensional 
wave is identically the direction of arrival of the incoming 
signal. 

(3) The estimation of the "phase" can be 
done relatively easily for the noisy trace function by mini- 
mum mean square error surface fitting. 

b. Linear Array 
The configuration of this array is shown in 
x Figure 1-7 
From Appendix A: 


- Time delay trace function 


AT.. Z S(i-j)sine (II-9) 
mI C 
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- Phase difference trace function 


k ھ0‎ . 
وره د‎ - AT Qy 7 EB c (II-10) 


where d - element spacing. 





Figure II-15 


Linear Array Configuration 


In Figures II-16 through II-20 are shown the three dimensional 


plots of the phase difference trace functions Ab. for a 


d] 
linear array of 32 elements (Equation II-10) of the same 
aperture as the circular array. The computer program is 
presented in Appendix B-2. 
Observations: 
۰ (1) The trace functions of a linear array are 
always a sampled two dimensional plane surface. 


(2) The "slope" of this plane is indicative 


of the direction of arrival of the incoming signal. 
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(3) The estimation of this slope can be done 
easily for the noisy trace function by minimum mean square 
error surface fitting. 

General Trace Functions 

For those arrays which do not have analytic trace 
functions the approach is slightly different. In this case 
in order to use equations (II-5) and (II-6) the value of the 
vectors m is needed either in cartesian (y) or polar 
(r,9) coordinates. Then each T; (q. II-2) is calculated 
separately as referenced to the origin of the coordinate 
system, and the differences are calculated and plotted. In 
this section two geometrical configurations were chosen for 
demonstration: 

a. A "Conformal array" which is composed of a half 
circle and two linear parts. The configuration is shown in 
Figure II-21. For this array, although it is built of well 
defined geometrical forms its trace function cannot be 
expressed in closed form. Figures II-22 through II-28 present 
the trace function for various look directions (bearings). 

It can be observed that it still preserves symmetry around 
the zero bearing and it is built up of pieces of the trace 
functions of circular and linear arrays. 

b. A "Random array" is illustrated which is composed 
pr 19 ت5‎ s randomly distributed within a circle of about 
9 meters. The configuration is shown in Figure II-29. The 
trace functions for various look directions (bearings) are 


Ehown in Figures II-30 through II-34. These trace functions 
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Figure II-21 


Conformal Array Configuration 


are totally random and only two obvious types of symmetries 
are preserved: 


E ۵9 د و‎ = pis 


- the whole trace function is antisymmetric for 
bearings of 180 degrees apart (this second 
characteristic is a result of the first one). 

The computer program which generates both the "conformal" 


and "random" trace can be found in Appendix B-3. 
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Random Array Configuration 
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Common characteristics of the non-analytic trace 
functions are: 

'" They do not have a distinct single parameter 
which is indicative of the bearing of the incoming signal; 
however, they still have a distinct form for each bearing. 

- As a result of (a) the minimum mean square error 
surface fitting still can be used; however, each trace func- 
tion must be fitted on a point-by-point basis. Some 
"educated" search procedures may be used. 

NOTE: 1. All the trace functions plotted in 
this part were phase difference trace 
functions. The time delay trace 
functions were not given because 
they are only scaled versions with 
the same functional form. 

2. Each intersection in the trace func- 


tion plots correspondsto a data point. 


E. DEGENERATE TRACE FUNCTIONS 

In some applications the task of generating and processing 
the full trace function may not be necessary or possible. For 
those situations a degenerate version of the trace functions 
can be applied. The degenerate. versions are derived by taking 
only some distinct pairs from all the possible pairs of ele- 


ments. 
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For example for a circular array a logical choice of 


pairs will be on opposite sides of diameters such that: 
a. N 
30 G t 7) mod N (II-11) 


Then, by substitution in (II-7) 


mee 2۱ e 
Bi 一 a COS (I 8) (II-12) 


and in (II-8) 


۵٩  ء‎ = AT..W, = —— cos (= - 8) (II-13) 
which is a simple cosine wave whose phase is identically the 
bearing of the incoming signal. It can be shown that the 
performance of the degenerate trace function is not as good 
as that of the full trace function but may be adequate for 


some applications. 


F. CHAPTER SUMMARY 

In this chapter the concept of the trace function was 
developed and demonstrated. From the few examples shown it 
can be concluded that the concept is well adapted to be used 
in a generalized spatial replica correlation scheme. Tn 
addition, for highly symmetrical arrays the amount of pro- 
cessing requires is significantly reduced. 

femelle cneseireular array 1S a natural one to be treated 


this way. In subsequent chapters the bearing estimator for 
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a circular array will be developed based on trace function 
concepts. In addition to the circular array the linear array 


Will be also treated mainly for comparison purposes. 
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III. MINIMUM MEAN SQUARE ERROR SURFACE FITTING AS 
APPLIED TO BEARING ESTIMATION BY TRACE FUNCTIONS 
A. INTRODUCTION 

In the previous chapter it was shown that the trace 
function can be represented by a three dimensional surface 
(function) whose shape is dependent on the bearing 0. For a 
given array the trace functions are known either analytically 
or in the form of a data matrix at any bearing 0. For the 
received data from an array of sensors the trace functions 
are estimated via a cross spectral estimation procedure, and 
are always in a data matrix form. In order to estimate the 
bearing by means of the trace function one has to "match" 
the estimated trace function to the known (stored or analytic) 
trace function. There are several ways to perform the above 
"match", some of them require the knowledge of the noise sta- 
tistics and others do not require it. As the noise statistics, 
especially the cross noise statistics, may be inaccessible, 
the methods which will be followed in this analysis will not 
require this knowledge. 

The best known method of this class is the minimum mean 
Square error (MMSE). In this method the sum of squared errors 
between the estimated (measured) value of the trace function 
End the "known" value of it at all coordinates "i,j" will be 
minimized with respect to the parameter "6" (the bearing). 

Ihe value of 8 which will give the minimum is the expected one. 

In order to give the reader a better intuitive feeling 


about the estimation problem Figures III-1 through III-10 
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show noisy trace functions one would get after the cross 
spectral estimation process. The above figures show trace 
functions for circular and linear arrays for various noise 


2 


variance levels from .001 Rad au dove J S6deg.) to 1.152 


Ed" (St.dev. 61.5deg.). It can be seen that even for the 
highest noise level shown one can still visually perceive 

the form of the trace function while for a single data point 
this noise may be overwhelming. The cause of this perception 
1s the visual integration over the surface of the trace 


function. The very same characteristic is used by an esti- 


mation processor. 


B. DERIVATION OF THE ESTIMATOR 
As it was shown the phase difference/time delay trace 


Function (46:4 and ER g) is dependent on indices "i,j" 
3 2 


8 
and parameter ۰ 


NOTES 0 will be used throughout the derivation but 


03006 


the analysis applies as well to ۵۲ د و‎ 6・ 
5 


The estimated trace function will be: 


fp (i,j,0) + v(i,j) Gem)‏ = ووفك 


where for 1s the trace function and v(i,j) is a random esti- 
mation error which will be approximately gaussian for long 
enough averaging time in the estimation process. As it was 
shown in [26] the error variance is not dependent on the 
estimated value. However v(i,j) has a covariance matrix of 


order (N 2xN2) which is complicated and its knowledge is not 
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assumed in the forthcoming derivation. Its’ effect on the 
estimator however will be:checked by simulation. 
Ihe general. 'weighted" least squared cost function is 


defined as: 


N-1 N-1 ^ ? 
E» >, Wi. (Aij = Abi; لم‎ (III-2) 
120 0 







The weights Des can be used to account for different 
quality of the measurements ۵9 و و‎ « and must have two main 


features: 


٣٠۶ ٢ر.‎ 
SCH Jr 


0 for izj 


W.. 
ij 


Minimizing the cost function J with respect to 0 will 
glve: 


oA $ 


ES (A 


)(-2 35 


N-1 N-1 
0 لا ل۷ دلو د‎ w. 
1220 0 
This is as far as one can go without specifying the trace 
Bunction explicitly. 
To solve the minimization problem two approaches are 
possible: 
- For "analytic" trace functions, to derive an explicit 


estimator. 
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- For "random" trace function to conduct an exhaustive 
or "educated" search (which is not efficient), or use 
some form of gradient search algorithm to achieve 
minimization. 

In this dissertation the first approach will be followed 
Memento the fact that the circular and linear arrays have 
"analytic" trace functions. 


| WW Circular Arrays 


The trace function for circular array is (Eq. II-8): 


1 0 cy. 1 RF š 1 
49 6 =p sin EEJ sin )اکسا‎ 


K sin pies] sin gi» (TIT-4) 


where K is a constant for a given trace function. 


Inserting the above equation in (III-2) and (III-3) gives: 


N-1 N-1 ^ T 1 2 
E y Wi; (401; - K sinly(i+3)-0] sinly(i-3)]) 
i-0 i-0 


(III-5) 
and 
N-1 N-1 - 
0 = ےت‎ 2K 之 を و و شاد وا‎ - K sin [i Ci*j)-0] sinfy(i-3)1) 


x sin[ そ (1 ュー] ) ] cos [7(itj)-6] (III-6) 
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0 S : Tes Aan ای‎ 
0 2K > W450i. Sin [Ti [ ( 1201 3) 6 ] 


- KW. .sin[(i+j)-01cos[T(i+j)-01sin [TCi-j)] 


1) 
N- N-1 
0 = X Wi ja9j 《sin 人 A ein) 4 
=0 3=0 
N-1 N-1 
" 27 2T 
- ECH (cosTi - cos; sinê 
1=0 j=0 
N- N-1 、 
K - 47. - 47. 5 - 27. 2T. 
-p 2 Wy Simp t sing) 2singpi cos q 
1=0 3=0 
-2cos i لکد زو‎ 06 
N-1 N-1 
K 4m 4m ZE: Es 
or WW,  (cosqpi + cos4rj- 2cosqpi cost] 
1-0 3=0 
+2sin ZT; sin^7j)sin20 (III-7) 


N 


In general the last two terms of the equation are non-zero; 
however, for some important and practical choices of Dr they are 
zero. 

Two cases for which the last two terms are zero are: 


-for Wij = 1 for all iZj 
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-for RS being equalin groups such that each pair which 
has the same value for cos LS (i-3)] will have the 
same W... 

1] 
The importance of the second case is due to the fact that 
۳۰۰ which are at equal distance in physical space will have 
the same weights. 
x For the cases when the last two terms of equation (III-7) 


are zero the estimator will have a particularly simple form: 


N 
A2. et aq 
۹ k: ت7‎ E v, bo; (sinri - singi) 
9. = tg "M o c > TTD 
N-1 N-1 
o W. 7 ¿(cosSpi - cost) 
i-0 j-0 


In the general case however, when the last two terms are 
not zero, the derivation is different. One can write the 


following equation from (III-7): 


K, cos9 - K, sinO - K. cos20 + Ky sin20 (III-9)‏ ے0 


ATK K, are expressions which can be easily 


۱۳۰ GRA 


calculated for a given estimated trace function. 


where K 


After some manipulations one obtains: 
v 


Daa K Y 1+tg“ - K 


1 tg0 


tgð y lttg'0 - K, (1-tg^0)*2K 


2 4 


(TII-10) 
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Which is an equation in one variable and can be written 


as 
0 = K, Y 1+x° - K2XY 1+x° - K, (1 -x°) 由 کا2‎ (III-11) 
and 
dc = tg x (III-12) 


The above equation can be solved by numerical methods. As 
mentioned before this procedure is required only if the weights 
are almost arbitrary. 

An additional feature of the estimator in (III-8) is that 
it does not require the knowledge of the constant K, which 


contains the frequency F, diameter D and speed of sound C. 


eS 2/۵ 
In fact it is important to realize that after estimating e, 
K can be estimated and as a result the speed of sound C, can 
be calculated when the frequency F and the diameter D are 
known. This feature is important for some applications such 
as seismic signal processing where there exists the problem 
of dispersive propagation in which the speed of sound is 
varying with frequency. The estimation of K will be done as 


Follows: 
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using the same cost function (Eq. III-5) but minimizing with 


respect to K when 6 is known: 


N-1 N-1 
A ES > 3 ;5 0 ach sinl= y(1+3)- ال‎ (i-j)]} 
120 jz=0 


ET] (III-13) 


End the estimator is: 


N-1 N-1 


ー > Wi ۳ ۰ sinl™(i+})-6] sin[gCi-j)] 





K (III-14) 
De ال‎ 
$ و‎ In). Gs 
Wij sin [で は 3 (۰۵ و[‎ 1 [gC Jl 
120 0 
knowing K: 
Peete (III-15) 


2. Linear Arrays 


Following the same development as used for the circular 


array, the trace function is (Eq. II-10): 





۳ 0 


13,0 C (i-j) sin 6 





= K, (i-j)sine Ciri=16) 





N-1 N-1 | 
2 F Been 2 
T = Hae [Ai د‎ - K, (i-j)sin®] (III-17) 
2۳-۱ 0 
and ` 
WEN N .-1 
BI ` ^ سر بج‎ E 
0 = wg -2K رز -1) [ 9126( -1) وک و و۵۵ اد و۷‎ 90 
120 3-0 
(III-18) 
and N-1 N-1 
> > inne) 
E el i=0 3-0 b 
di Sin N-T - (III-19) 


120 3-0 


3. The Degenerate Trace Function Estimator 





The degenerate trace function for circular arrays is 


EG IL-13): 
i ES DE 27 
Biz م‎ = a cos (mi 8) (III-20) 
where it was assumed 


jJ = (i+ N/2) mod N CT iene? 1, ) 
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Inserting this equation to the estimator (Eq. III-8) 
and assuming all re equal to unity: 
KEN 
٨ ۳ 
Ad; ا‎ a 
۸1 6 کت ےڈ‎ — ——- TI) 


نه وه | > 


which is the phase of the first harmonic of the discrete 


fourier transform of Ad. . 


C. CHAPTER SUMMARY 

In this chapter the MMSE estimator of the bearing 
was developed for circular and linear arrays based on the trace 
function principle. Two general remarks should be mentioned: 

l. The . MMSE estimator is equivalent to the Maximum 
Likelihood estimator when the noise is gaussian and indepen- 
dent among the data points. This case can be considered the 
"best" case and it will be used to develop the  Cramer-Rao 
Lower Bound for the estimator. This CRLB . will serve as 
a benchmark to the performance of the estimator which was 
derived by simulations. 

2. If all the data points in a trace function are esti- 
mated from the same time frame/s the two halves of the trace 


^ 


function are totally dependent, because ed 70 6 





from the same 'measurement". In addition the terms ۵4 are 


fnot contributing to the bearing estimation process. 
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As a result of these the number of contributing terms in 


mine trace function data matrix is: 
(N 2 - N)/2. 
The summation limits in the estimators must be changed 


accordingly: 


Equation (III-8) becomes: 








N-2 N-1 
y Nou. A 
2 Ree Ech 
A d -1 =0 J=++1 
0, = tg N= TT D'OR 
2 u: A. 
ات تا‎ 9 
| 1-0 ] = ユ +1 
Equation (III-19) becomes 
N-2 N-1 
) : Wi 5615 (1-3) 
B m1 0 1Z1+1 
9 = sin ° FET TNT )111 اه‎ 
Wis (1-3) 
i-0 J=1+1 
x Equation (III-22) becomes 
11/ 2-1 
2 2 0 
AQ. Ee 
3 ۲ -l 120 = 
96 p. = tg ` وو‎ 2 (III-25) 
0 ۱ 2 2T. 
Ao. cosi 
i-0 
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TUIS? - 






However, if the measurements for EE an PIS are done 
from independent data blocks the equations should remain in 
the original form. 

The exact structure of the weights Wi ュ was not treated 
here and it will not be pursued in this work, however, it 
should be emphasized that this is a subject for future work. 
Especially it is recommended to check the derivation of the 
weights with respect to their connection to the coherence of 
the two channels participating in each pair. 

In the next chapter the performance of the estimators 
which have been derived in this chapter will be checked by 
simulation to demonstrate the effect of changes in various 


significant parameters on the estimator performance. 
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IV. ESTIMATOR PERFORMANCE 


A. INTRODUCTION 

In the previous chapter bearing estimators for circular 
and linear arrays were developed using trace function prin- 
۰۳۱ 

In this chapter the performance of the estimators will 
be checked by simulations. The two main measures of perfor- 
mance are: 

- Bias 

- "Beamwidth" which here is defined as twice the standard 

deviation of the bearing estimate (20) 

In addition to the above a phenomenon called threshold 
effect was observed which is expected to be present in most 
nonlinear phase estimating processors [38, p. 287]. This 
x threshold was also checked for various parameters, and it was 
defined arbitrarily as the input phase difference variance 
level for which the output "beamwidth" is less than five 
degrees. 

The input to the estimator was the "estimated" phase 


difference .trace function (Eq. III-1): 


Ad,;, = (1,350) + v(i,j) (IV-1) 
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Two noise models were used for v(i,j): 

- independent gaussian noise from pair to pair 

ー correlated gaussian noise from pair to pair 
The model of the correlated noise is defined in Appendix C 
along with the Eu E method and computer program. The 
intent of including the correlated noise was to give a 
feeling about its effects on the performance and not so much 
on its accurate modeling. 

The simulation runs were performed mainly for circular 
arrays, however some performance measures were checked for 
linear arrays as well as for comparison purposes. 

ES een par meten ehange 48 computer runs were done in 
order to get reasonable statistics. The mean and the variance 
of the results of the 48 runs were calculated and they served 


as the estimated bearing and bearing variance. 


B. CRAMER-RAO LOWER BOUND FOR THE IDEALIZED ESTIMATOR - 

As a benchmark for the performance of the estimators the 
Cramer-Rao Lower Bound (CRLB) was used. It was calculated 
for the maximum likelihood estimator which is equivalent to 
the MMSE estimator for uncorrelated gaussian noise 
[23, p. 198]. This case can be considered as "best case" and 
it can be used as a lower bound for any realizable estimator. 
The CRLB for a maximum likelihood estimation of a parameter 


IS eh. p. 275. BW p. 335). 


2 x No 


9s(t,a) 2 
00 


(IV-2) 


T 
27 dt 
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where No - the input noise level (variance) 
s(t,a)- the "clean" signal 
o ` - the parameter to be estimated 
Adaptation of the above to the present case will be: 


2 
- No = o 22. - variance of the estimated phase difference 


Ad. - 


1] 
- S(t,a) is the discrete trace function f. (i,j,0) 
- The integration becomes double summation 


and for a circular array: 
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The above bounds were plotted along with the simulation 
results for comparison purposes and, as it will be shown the 


simulation results achieve this bound asymptotically. 


en ) 5ض‎ 
1. Circular Array 
For this performance test the simulation was run for 
4, 8 and 16 element 5 meter diameter circular arrays at 100 Hz 
and target bearing 30°. For each element number change, the 


input phase difference variance was changed between ۱ , 0 * 


5 2 


(120°St.Dev.) down to 6x10 ? Rad“ (.45°St.Dev.) in 25 steps. 


The bias was calculated for each run by the equation: 


B = 6 - ۵ 
where B - bias 
0 - estimated bearing 
0 - peal bearing - 0 


The results are plotted in Figure IV-1. 

The "lock-on" range was marked on the plots and although 
the bias was not the criterion for the "lock-on" range defini- 
tion, it is clearly shown that the bias fluctuation becomes 
very small within the "lock-on" range. 

Observations: 

a. Within the "lock-on" range the bias is very small 
and is bipolar, a fact which suggests that the apparent bias 


is due to an insufficient number of runs for: ensemble averaging 
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and that the estimator itself is really unbiased. 

b. Even if "a" is not true the estimator is at least 
asymptotically unbiased. 

c. The correlated measurement noise has no apparent 
effect on the bias compared to the uncorrelated noise. 

d. As the number of elements increases the fluctuations 
in the bias are smaller. 

2. Linear Array 

For this performance test the simulation was run for 
a linear array of 8 elements with 1.91 meter element-spacing 
at 100Hz (taking a 5 meter diameter circular array with the 
same number of elements and opening it up). However, this 
time the bearing was varied as a parameter because unlike the 
circular array the linear array is not symmetric in the x-y 
plane. The range of the input was the same as for circular 
array. Figure IV-2 shows the resultant plots. 

In addition to the observations of the circular array, it 
can be seen that the bias fluctuations increase toward end- 
fire but even at those bearings the bias is asymptotically zero. 

9. Conclusions With Respect To Bias 

a. It can be assumed that the MMSE . estimators 
developed in this work are unbiased or at least asymptotically 
unbiased. 

b. The lock-on range is a good measure for the region 
where the bias fluctuations settle down. 

c. Increasing the number of elements reduces the 


bias fluctuations. 
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D. BEAMWIDTH 

The beamwidth, which was defined as twice the standard 
deviation of the bearing estimates (208), is by far the most 
important performance measure. As a result it was most 
extensively checked both by varying the various parameters 
which effect the performance and comparing to the CRLB. In 
this section most of the simulation runs were performed for 
both correlated and uncorrelated "measurement" noise along 
with the CRLB. 

The model used for the correlated measurement noise is 
presented in Appendix C. This case of correlated noise is 
the actual one when the phase difference estimation for all 
pairs is done simultaneously (from the same time frame). 
However, if the estimation will be done from different time 
frames for various groups of pairs the measurement noise can 
be accounted as uncorrelated and then the uncorrelated mea- 
surement noise case will apply. The effects on the perfor- 
mance of both cases is shown throughout this section. 

l. Circular Arrays 


a. Beamwidth Versus Phase Difference Estimation Noise 
Variance 


Figures IV-3 through IV-6 show the results of the 
Simulation runs for this performance measure for 4, 8, 12 and 
16 element circular arrays of equal diameter of 5 meters at 
100۳2۰ For the 4 and 8 element arrays both correlated and 
uncorrelated measurement noise is given. For all the plots 
the CRLB is shown along with the simulation results asa 


benchmark to performance. 
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Observations: 


De 


In all plots the results and the CRLB exhibit a 
clear threshold effect. 

The correlated noise assumption (which is the 
realistic one) does not have a major effect on the 
performance. 

Under both correlated and the uncorrelated measure- 
ment noise assumption the estimator achieves 
asymptotically the CRLB at low level input noise 
variance, thus the estimator is asymptotically 
efficient [38, p. 276] for decreasing input noise 
level. 


Beamwidth Versus Numbers Of Elements At Constant 
Diameter 


Figure IV-7 shows the dependence of the beamwidth 


on the number of the elements for a circular array of 5 meter 


diameter at 100Hz. The plots were done for three representa- 


tive input variance levels. Only the uncorrelated measure- 


ments along with the CRLB are shown which are adequate to 


F present 


the basic behavior. 


Observations: 


- A saturation effect is clear for all three variance 


levels shown; at each noise level at a given diameter 
and frequency there is a limit on the number of 
elements required, increasing the number of elements 
above this "saturation number" will improve perfor- 
mance only marginally. For example for the lowest 


3 


curves (o 2 =2.25x10 Rad“) this number is 
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approximately 16 elements, for higher noise levels 
it will be higher. 

- Again the simulated estimator performance approaches 
the CRLB asymptotically for large number of elements. 


c. Beamwidth Versus Diameter For A Constant Number 
Of Elements 


Figure IV-8 shows the dependence of the beamwidth 
on the diameter of a circular array of 8 elements at 100Hz. 
The plots were done for two representative input variance 
levels. The correlated and uncorrelated measurements along 
with CRLB are shown. 

Observations: 

- A saturation effect is clear in this case too; 
increasing the diameter of the array above certain 
limits will improve performance (smaller beamwidth) 
only marginally. This limit appears to be approxi- 
mately at 15 meter diameter which is one wavelength 
at 0 2 ۰ 

- The behavior of the simulation resultant curves is 
Similar to the CRLB and they approach the CRLB 
asymptotically. 

2. Linear Arrays 
The analysis of performance for the linear array is 
Beene lcms extensive than forfthe circular array and its aim 
is two fold: 

- It is the most standard geometrical configuration 

treated in literature and a reader of this work 


can readily compare it with other results. 
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- It was compared to the circular array in some 
aspects of its performance. 

In this section only beamwidth versus phase difference 
estimation noise variance plots are presented. Figures IV-9 
and IV-10 show the resultant curves for a linear array of 
5 meter length (same aperture as the circular array) with 
4 and 8 elements at 100Hz. The plots are for both correlated 
and uncorrelated measurements along with the CRLB. 

Figure IV-11 and IV-12 show the results for a linear 
array of 8 elements and 1.91 meter interelement spacing. 

This configuration is equivalent to an "opened up" version 

of a 5 meter diameter circular array preserving the inter- 
element distance and number of elements. Figures IV-ll and 
IV-12 differ in the look direction which causes a different 
performance, as expected, because the asymmetry of the linear 
array in the x-y plane. 

Observations: 

- The linear array estimates present a similar be- 
havior as the circular array with respect to their 
performance relative to the CRLB - they are 
asymptotically efficient. 

- A linear array with the same aperture and number 
of elements as a circular array has a slightly 
worse performance, the reason being the trend to 
concentration of elements at the ends for the 
circular array, a characteristic which improves 


bearing estimation. 
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~ A linear array with the same number of elements 
and linear-element distance has however a much 
better performance than the equivalent circular 
array because of the much bigger aperture. 

- For a linear array the beamwidth increases toward 
endfire, a fact which is expected from the equation 
of CRLB (Eq. IV-4) where the dependence on 1/cos “8 
is observed. The simulation results show the same 


behavior. 
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NN. THARESHOLD EFFECT AND ITS IMPORTANCE 

As mentioned in the introduction to this chapter a 
threshold effect which can be expected in this kind of non- 
linear estimator was observed. In this work the threshold 
was defined as "lock-on point" and it is the value of the 
input variance for which the output beamwidth becomes less 
than five degrees. The number of five degrees was chosen 
arbitrarily and it can be any number in the vicinity of the 
knee of the performance curve (see for example Figures IV-3 
through IV-6 for circular arrays and Figure IV-9 through 
IV-12 for linear arrays). 

The range of input variances for which the beamwidth is 
less than the threshold value, in this case 5 degrees, is 
defined as "lock-on range". 

The importance of this issue is due to the following: 

- The lock-on point can be translated to a threshold on 
the coherence at a specific frequency and only those 
phase difference measurements for which this threshold 
1s exceeded are processed by the trace function estima- 
tor. (This connection will be further explained in 


Chapter V.) This thresholding may ensure that only 
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those frequencies which will guarantee a certain beam- 
width will be processed. It should be mentioned that 
the threshold is frequency dependent exactly in the 
same manner as the diameter dependence which will be 
shown later in this section. 

- As it can be seen from Figures IV-1 and IV-2 operating 
within the "lock-on range" ensures also that within 
the beam width there will be less variability of the 
resultant bearing estimate. 

The dependence of the lock-on point for a circular array 

on number of elements and diameter is shown in Figures IV-13 
and IV-14 respectively. The plots are for uncorrelated 
measurements and for the CRLB. 

The frequency dependence of the lock-on point is exactly 
the same as for diameter and is not shown separately (both 
have an "xn effect). 

The values of SNR shown on the plots are obtained from 
variance values related back to the system input and their 
aim is to give a feeling of the change in the noise tolerance 
of the system as the diameter (or frequency) and number of 
elements are changed. 

Although the simulation results do not achieve CRLB 
within the range of values presented, nevertheless the trend 
shown suggest that it is eventually asymptotically achieved 


for higher number of element and larger diameter. 
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F. PERFORMANCE SUMMARY 

1. Within the "lock-on range" the estimator is practically 
unbiased. 

2. With respect to the beamwidth performance the estima- 
tor is asymptotically efficient i.e. it asymptotically achieves 
the CRLB when the input variance decreases, or the number of 
elements, diameter, or frequency increases. 

3. A threshold effect, defined here as "lock-on" was 
observed and it enables one to tailor the desired system per- 
formance. 

4. There is a limit on the number of elements at a fixed 
diameter (or on the diameter at a fixed number of elements) 
beyond which the beamwidth will only improve marginally with 
the increase of number of elements or array diameter respec- 
tively. 

5. The trace function estimator applies to circular 
arrays as well as to linear arrays, however for the same 
aperture and number of elements the circular array performs 
slightly better. 

In the next chapter on system considerations a system 
which uses the estimators developed here will be suggested. 
The problems involved in implementing such a system will be 


analyzed. 
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V. SYSTEM CONSIDERATIONS 


A. INTRODUCTION 

The motivation behind this work has been to find a solu- 
tion to the practical problem of improved bearing estimation 
using a small array. The trace function concept offers a 
way to solve the problem but is not the solution itself. In 
order to solve the problem giving a meaningful and practical 
result the "pieces" must be connected together into a system 
concept which considers the signal from the sea to the display. 
In this chapter such a possible system will be suggested and 
outlined. It is not the intent to design a system but rather 
to highlight those points which are to be considered in a 
system using the trace function concept. 

The main task would be to process the received data at 
the sensors (hydrophones) to bring it to a trace function 
format. 

Once the data is in phase difference and/or time delay 
trace function format the estimation process developed in 
previous chapters can be applied and the bearing estimated 
for each frequency resolution cell of interest or for the 
wide band signal. 

In the process of performing this task there is a major 
difficulty to be overcome; any processing scheme used must 
preserve the relative phase between the signals received the 
Same as they were at the receiving sensors in order to result 


in a meaningful trace function. 
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Two process steps may be identified which pose this 
11 17۷712 

- The decorrelation (prewhitening) both in space and time. 

- The phase unwrapping which is needed because of the 

inherent 27 ambiguity present in any phase measuring 
system. 
Those and other points will be reviewed and analyzed briefly 
in this chapter. Some suggestions will be made to overcome 
the mentioned difficulties. 

In Figure V-l a conceptual system functional block dia- 
gram is presented. It applies the trace function principles 
as discussed in previous chapters. In the forthcoming sections 
the functions and problems associated with various blocks will 


be treated. 


B. SENSOR ARRAY 

The sensor which will usually be a hydrophone is the first 
part of the system interacting with the signal. The major 
characteristics which have to be considered at the sensor 
array are: 

- Accurate and known geometrical configuration. 

- Equal phase response of the sensors across the frequency 

band of interest. 

If the phase response characteristic cannot be satisfied 
because of practical problems, at least an accurate response 
chart at all frequencies of interest must be obtained and 
stored in a computer memory. This data can be used in turn 


to compensate for the different response characteristics of 
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the various sensors when the phase difference calculations 
are performed for each pair of sensors. This way the rela- 
tive phase differences at the input to the system can be 
preserved. Any other change in response of the sensors will 
be treated by the system as additional noise on the phase 


difference measurements. 


C. PREPROCESSTNG 

In this part of the system preamplification, bandlimiting 
and A to D conversion takes place. The only requirements are: 

- Zero phase or at least linear phase response at all 

frequencies of interest. 

- Equal phase response for all channels. 

After the A to D conversion digital methods will be used 
and the preservation of equal phase response becomes much 


simpler. 


D. SPATIAL DECORRELATION (WHITENING) 

The processing of the signals received by an array of 
sensors in a noisy background is based upon the relative 
coherency of the directional signal compared to the relative 
incoherency of the background (ambient) noise. When this 
difference is distinctive (i.e. the noise is almost uncor- 
related from sensor to sensor) the processing is relatively 
easy and not very much dependent on the noise model. This 
characteristic is achieved when the separation between sensors 
15 in the order of a half-wavelength for the frequencies of 


interest. The task is very much complicated where the separation 
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is small, because the high correlation of the background noise 
narrows the difference between the signal and noise in that 
domain. 

In the problem defined in this dissertation the case is 
emer. arrays where the whole array aperture is much smaller 
than a wavelength at the lowest frequencies of interest and as 
such the noise will be highly correlated. 

The conventional way to decorrelate the noise among the 
sensors is by multiplication with the inverse noise covariance 
matrix or the inverse signal and noise covariance matrix. 

This operation is prescribed by every "optimal array processor" 
(2, 5, 39) and it can be performed based on a priori know- 
ledge," estimate and plug' or adaptive calculation of the in- 
verse matrix. 

This operation is inherentely not phase preserving and 
this fact will be shown in the forthcoming simple example: 

Assume two sensors with the signal part of the input 
being: 
سے‎ 


x = = (V-1) 
oi fwttAd) 


= 


27 


clearly the phase difference is Ab. 


A typical covariance matrix has the form: 


1 p 
و‎ )۷-2( 
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where 
2 
= ーーー ——— — 
haba A) 
The prewhitening process is defined as Q^*x as it can be seen 
from: 
: EL __ = = ER = یا ا‎ 
E {0 "x xt Q 3] = Q 5 E(x x^) Q E AO OS => TE 
and 
1 -a ۱ ۱ 
E 1 elut 1 ال یج _تالالی‎ ( 
Q fx = —T EE 
ch -4 1] تا‎ | .lwt, رو‎ 
(V-4) 
iot S lut : 
A _ e " 140 ~ Te LAG E 
x, = — (1 ae ) هن‎ ee (e a) (V-5) 


After some manipulations the phase difference between 


2 l 
2 ۳ -1 5 ٨ -1 sind6(o2-1+ A-02) 
۵ x = tg EE Kë EE D Di سے‎ 
21 cosAd(l+ta”)-2a (cosA6-0) (1- A -p2) 
(V-6) 


which clearly cannot be equal to Ad for any value of p except 
p>0 which the trivial uncorrelated case. 
In order to preserve phase from the total system point of 


view one must correct for the above phase distortion. 
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A new method called "phase preserving spatial decorrela- 
tion (or prewhitenning)" by the author is proposed and out- 
lined here. The notion is that the phase distortion caused 


by the decorralation process can be calculated and compen- 


sated for, after the cross spectral estimation process, pro- 


vided that one knows exactly the inverse covariance matrix. 
To demonstrate the above, the previous example is con- 
tinued: 


the original phase difference was ۵۹ 5 


9 


the resultant phase difference after decorrelation was 


E : 2 
MM- ~ = ایم‎ 5 1۳۵4 )1-2 ) 


"uto cosA ゅ (1+a^)-2a 


The difference between the two is the compensation term: 


-1 2asinAd(1l-acosAb) 
l-a 052۵4 -2 4 


۸۵-۵۵ = tg (V-7) 


The above compensation term is impractical to calculate be- 
cause it is dependent on the real phase difference which is 


unknown. As a result no simple compensation can be used. 


Equation (V-6) can be solved numerically however for Ab, 5 
۳ JE 
once A ゅ て 9 is measured provided p is known. 


2 


As shown by this simple example this approach is very 


2 


complicated even for simple arrays (2 elements) because it 


has to be calculated for each frequency bin and each sensor 
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pair. For practical number of elements in an array this kind 
of operation becomes prohibitive. To solve the problem a 
different approach has to be found. 

As the ambient noise is a wideband process its cross 


covariance function has a limited time width 
رہ‎ NEBOT > ie (V-8) 


i.e. it can be made arbitrarily small if long enough delay 
between the elements is applied. 

Figure V-2 shows such a covariance function for a band- 
width of 50 to 250 Hz a power spectrum of -6dB/octave for 
2.5 and 5 meter distance between sensors. 

The graph clearly shows the decay of the covariance "p" 
as function of delay time "t". This characteristic can be 
used for the decorrelation process. One can introduce 
successive and known delays between the elements of an array 
such that the smallest delay is enough to decorrelate the 
noise down to a desired level for the closest pair in the 


array. 


After this process the signal is: 


(E) 
0 
RC) = (V-9) 


Xy] (t-(N-1)4) 
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Spatial Covariance 
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Figure V-2 
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where A is the incremental time delay introduced between the 
elements. 

It is obvious that to compensate for this process the 
only operation to be done is to re-add at each frequency 
We a phase compensation term to the calculated phase dif- 
ference. This term will be for the pair "i,j": 


Ad = ACi-j)w, (V-10) 


BJ comp. 
This method will not affect narrow band signals but may also 
decorrelate partially the wide band signal. 

This effect can be minimized if A is chosen such that is 
large enough to decorrelate the noise but (N-1)A is smaller 
than the decorrelation time of the wide band Signal. 

In order to perform this operation one must have the 
moise and signal correlation function which can be obtained 
either from apriori knowledge or from estimations made in 
situe 

The above method should be viewed as a suggestion and not 


a firm solution mainly because it is based on theory only 


without simulation or measurement back-up. 


E. REVIEW OF SPECTRAL ESTIMATION 
This section will summarize the spectral estimation tools 
needed for subsequent estimation of time delay and phase 


difference. 
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As it will be shown later if one has infinitely long data 
record he can get both stable estimates and infinitely good 
resolution of the power spectra, and phase spectra of a 
stationary random process. The main factor which limits 
achieving the desired performance is the limited length of 
data available. The reasons are many, but one of the most 
important ones is the time limit on the stationarity of the 
environment. In the case of data in the ocean, a few minutes 
can be considered as a good measure for stationary environment. 

The analysis which will follow reproduces the main results 
in Ref. [25, 26, 27] which deal with the limited time record 
auto and cross spectral estimations. 

The two basic parameters which control the performance 
of spectral estimation are the available record length, T, 
over which the sample of the random process is assumed 
stationary and the desired frequency resolution, B, of the 
spectral analysis. Large values of BT give good performance 
Dut small values are often encountered because of short T 
available or requirement for fine resolution B. 

In order to take advantage as much as possible of the 
available data a method was developed [40] to obtain good 
reduction of variance by overlapped processing of the data. 
This method is less efficient than the same number of segments 
non-overlapped but it is much better than using the available 
data in non-overlapped form. 

Several questions arise when one attacks the problem: 


what percentage of overlap to use?, which windows?, how big 
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the FFTs have to be? Those questions and others were answered 
in Ref. [25, 26] and the results follow. 

The stationary random process x and y are observed for a 
time T seconds, thus x(t), y(t) are available for O<t<T. One 
Wishes to estimate the cross spectrum Gy Sf) with a resolution 
of B Hz over a bandwidth of Q Hz. 

The method of ar و‎ the available data is shown in 


Figure V-3. 


w Data Window 


AE 





O S L S-L (P-1)S (P-1)S*L 


Figure V-3 
Segmenting Data With Overlapp 


The connection between the various terms in Figure V-3 is: 
(P-1)S*L « T (V-11) 
where P - number of segments 
S - shift of each adjacent window 
L - data window length 


T - total available data length (time) 
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1. Estimation Procedure 
The estimation of the auto and cross spectrum is done 
in the following way: 


a. DFT the segmented and windowed data x and y: 


SE 
X CE) = X.(3) wi C) e" i?7jn/L ( VS 12) 
] =0 
Esch 1 د‎ 
ei, YG) wG) و‎ (V-13) 
(0 
where 


X (f ), Y (f ) - Fourier transformed data at frequency 
k n K m nF 
E 5 : 
f ہت‎ of segment k 


- time index (discrete) 0,1,2...L-1 


k - segment index 14 17.2 
n - frequency index 0 +2 
F. - sampling frequency 


NOTE: The index j begins from zero at each segment k. 


b. Estimate of auto power spectra 


l P 
E Ix. CE 9]? (V-14) 
mx n E ۳ 1 
s ^ 
=] 
E 
ED “EFT 2 IY C£.) | (V-15) 
& KEI 
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c. Estimate of cross spectra 





P ^ 
۸ E > # 2 TER CE) 
k=] 


Tm[G (r) 1] N 
n n | = ۵۵ (f ) (V-17) 


Phase Ph, |٢ ee) 
sa Re[G (f )] 
xy n 


RENNES A ce) = Rele (f )])?+(Im[G (EIN? (V-18) 
xy n xy n Xy n 


d. Estimate of magnitude squared coherence 


^^ ; 2 
A | IG (en 
د )نا‎ — (V-19) 
Ee E) 
NOTE: The magnitude squared coherence function whose 
role will be discussed later, is given here because 
of its connection to the other estimation processes. 


Be DE The Estimates 


The spectral estimates presented are random variables. 


Their means and variances will be stated under the assumptions 


that x and y are Gaussian random processes and the frequency 


resolution of the spectral window ۱۷| 7 is narrower than the 


finest detail in spectrum G. 


where 


L-1 
we) = 2 w(3) e 1273n/L (V-20) 


(0 
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a. Means 


^ e 2 
E (6, (£f )) Sa f Ilwv) | “ay (V-21) 


11 


E (6, (f) 


2 
i تہ‎ eo |WCv) | “dv (V-22) 


One can assume without loss of generality that the windows 
are normalized i.e. f |w | 2 = 1 then: 


^ 


E (G (f )] 
XX n XX 


00 
02 


(f ) (V-23) 
n 


22 
一 
Fh 


E 81پ‎ 29 (V-24) 


G 
xy 


Thus the estimators are unbiased 


E سا‎ (f. )) = Aby fn? (V-25) 
E fA, fn?) ? IET ا‎ (V-26) 
E v, 1^) = |y|* + s(1-|y る ) لپ‎ (V-27) 


(acer P28) 


The estimate is biased by the second term of the 
equation. A note about the above bias is given below. 

It is calculated for nonoverlapping segments. However, 
when overlap occurs this bias is substantially reduced to a 


value of about 50% of the above figure at 50% overlap. 
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b. Equivalent.number of degrees of freedom 


The convenient way to describe the stability of 


any positive estimate is its equivalent number of degrees of 


freedom (EDF) [18]. The definition of this stability constant 


"s: 


EDF = K = 2Laverage]* = CD 


(V-28) 
variance VE GCC 


Basically using tables as in Ref. [18] one can 


Berimate his confidence in the measurements. 


For example for K=500 90$ of the measurements will 


be in the range of 0.92 to 1.081 of the average.  Maximizing 
those K's will be one of the important criteria for the over- 


lap and window choice. 


c. Variances 
(P-1) 5 
2 : 2 1 k 2 
Var {G (f )) ミ G (f) i > 2 7 (V-29) 
k=-(P-1) 
q 
where 
$ (て ) = f w(t)w (t-T)dt (V-30) 
(BE) 
3 7 2 1 
Eu e ریت‎ » abel) po as? (V-31) 
K+ (P-1) 


At this point we can define EDF (k) for both auto and cross 


spectra [25, 26]. 


1.289 





28 


Ka = p 1 for auto spectra (V-32) 
$ (kS) 
Ca 
P $ (0) 
k=-(P-1) 
and 
5 2 
E- .ا‎ Ko for cross spectra (V-33) 
۹ 2 
= + = 
Var tA, (£3 Bg FG, (FEI E? ]/K, (V-34) 
の 
Varía9_ (f )}= —“L2— (V-35) 
P n |y... C£ MK 
xy n A 
Those are based on the assumption that 
و‎ 1 << 1-۰ It should be noted that neither of the two 


cross variances depend on the actual phase and the phase 
variance is not dependent on actual data values at all, but 
only on the coherence between the two sequences and the 
processing scheme. 


Var (|, (£, |^) 3 اح ان ا‎ for y, #0 


(V-36) 
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3. Design Tradeoffs 

As one can see from all the above statistics, k and P 
are playing major roles in the variance reduction. 

In order to improve performance one has to maximize 
k and P (which are connected) for a given amount of available 
data. 

In Ref. [25] this subject was thoroughly investigated 
for various window functions. 

Following is a summary of important results: 

a. The EDF(k) has a maximum as related to P fora 
given amount of data T. 

b. This maximum was proven to be almost independent 


of the window employed for large BT products 


max EDF * 3(BT-1) (V-37) 


c. The required fractional overlap which is defined 
as 


ae = EE D (V-38) 


is significantly increased for max EDF compared to 0.99 max 
EDF soit is not recommended to design for max EDF, but 1% or 
2% less. 

As an example using cosine window (Hanning) and 
BT=64 for max EDF P=146 segments are needed but for 0.99 
max EDF only 112 segments are needed. Clearly the additional 


30% of processing is not worth the one percent increase in EDF. 
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The optimal FO comes out to be in general not a 
convenient number to use, e.g..0.61 for cosine window. One 
will probably choose the closest convenient value like 
FO=0.5 which will give for cosine window EDF=92% of max EDF 
which is a reasonable figure. 

d. As the max EDF is practically independent of the 
window employed no tradeoffs had to be made to realize better 
side lobe response of the more sophisticated windows. This 


issue is effecting the size of the particular FFT's to be 


performed. 
E 
No ニ Bit (V-39) 
where 
N. = number of samples in one FFT 
C. = half power bandwidth constant of a window 
At = sampling interval 


Clearly for a "better" window C. is bigger in 
general and N. will be bigger. Thus the side lobe performance 
for a given B can be improved by paying with bigger FFT's and 
not affecting stability of the estimates. 

4, Effect Of Pure Tones (Lines) In The Spectrum 
The results shown so far were true under the assumption 
that the B of the spectral window IgE )| * is narrower than 
the finest detail of the spectrum. For tonal signals this 
Will not be the case and a change in the statistics will occur 
but the processing method will not change. In 
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Ref. [27] those changes were made and the following statistics 


derived: 
a. Mean 
EtG(f_)} š CIE) a8 + G(f ے‎ He (V-40) 
where 
G(f میدن(‎ = The previous wide band power spectrum 
GE = A^wC, -£ )|2 (V-41) 
n line 1 9 
A - line amplitude 
چ‎ line frequency 
2 u 2 
E (6, (Cf 2) = د‎ En uiae * AV IW CES 7E, | (V-42) 


D. "Variamee 


Var (G(f )) = Var (G(£, 2), 1a, *26 (f D, a Slt ) 


n narrow 
P-1 
x i ユー) 。 (ks )。12Tk(f -f。)S (V-43) 
P P W 
7 


from which a bigger absolute value of the variance can be 
"expected when a tonal is present. However, if the EDF (K) 


is calculated for Gtf_) >> ( Line which is the case in 


ide 
general when lines are present, one can see that this relative 


measure is increased which means a more stable estimation. 


1:38 





SEN 0005 7 Hue وی‎ 


r í و‎ 
EDF = K b n narrow 


tonal P-1 
Ge(f ) . > bla giel د‎ 


n wide 


k=1-P 
(V-44) 


this Ke onal 1s bigger than the one for wide band 


alone by a factor of approximately 


gG(f ) /G(f ) 201 (V-45) 


n narrow wide 


All the other equations for cross spectral estimation statistics 


۲ 5 111 hold with the new K 1 
tonal 


It is worthwhile to note the |y P will be much 


Xy 
bigger in this case (very close to unity) and its bias and variance 


will be very small because of the increased "S/N". 


F. PHASE DIFFERENCE ESTIMATION AND THE UNWRAPPING PROBLEM 


The basic phase difference estimation is done entirely by 


the cross spectral estimation process as reviewed in the pre- 


| 


۱ 


vious section. However there are some additional steps to be 


performed in order to construct the phase difference trace 
function or estimate the wide band time delay [32, 36]: 
- The frequency bins for which further processing will be 
done must be detected. 
- The phase difference estimated value must be recompen- 


sated for the decorrelation process as stated in Eq. (V-10). 
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| 


- The resultant phase difference has to be unwrapped 
because it is measured only in the range +n while it 
can have any +n2q ambiguity. 

1. Detection Problem 

The detection of the frequency bins which are candi- 
dates for further processing is important for the following 
reasons: 

- To reduce the false alarms in the system. 

- To reduce the computation load in subsequent stages 

of processing. 

Probably the most representative parameter of the 
phase difference integrity at a specific frequency is the 
に 


magnitude squared coherence |y This fact can be seen 


xy 
from Eq. (V-35). While KA is a fixed system parameter 

1 (7 is the data dependent term, or more precisely 
dependent on the ratio of the coherent to uncoherent signals 
between the two sensors considered. Moreover it can be 
related to the input SNR under certain assumptions of equal 


and uncorrelated noise between the sensors and linear pro- 


cessing up to the estimation stage, by the equation [6]: 


رد تج 
gu mir )۷-6(‏ 
Yxy n‏ 
As such the magnitude squared coherence PE can be a good‏ 
measure for detection [7]. There are two different processing‏ 


paths, each with slightly different requirements from the 


detection process: 
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a. Time Delay Estimation Path 

This path will require probably a lower threshold 
on the coherence because it is estimating the time delay from 
a large number of frequency bins, as it will be shown in the 
next section, and it further reduces the noise level, before 
the trace function processing. 

b. Phase Difference Path 

This path will have a higher threshold because it 
is intended to process mainly line or narrow band components 
which if they are present exhibit a fairly high coherence. 

To preserve the integrity of the trace function 
which includes all the pairs in the array it is not enough to 
detect a potential frequency bin on a "per pair" basis. This 
operation has to be done across the whole array (all the pairs). 
Two possible ways to perform this task are: 

- At each frequency bin each pair is thresholded for 
its magnitude squared coherence. Then the number 
of pairs which pass the threshold are counted and 
if this figure exceeds a predetermined number this 
frequency bin is qualified for further processing. 

- At each frequency bin the sum of the magnitude square 
coherence of all the pairs is calculated and if this 
sum exceeds a threshold this frequency bin is quali- 
fied. 

It should be noted that this detection process based 
on the coherence is a data adaptive process. It qualifies the 
various frequency bins based on its estimated magnitude square 


coherence in real time. 
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2. Compensation 

This operation is straight forward as explained in 
Section V-D. 

One has only to add the compensation term Eq. (V-10) 
to each "qualified" phase difference measurement according to 
its frequency and the participating sensors in the pair. 

3. Unwrapping Problem. 

This problem is well explained in Ref. [32]. There 
are several solutions treated in the literature [32, 36] how- 
ever in the case of the circular array or similar two dimen- 
Slonal arrays the problem can be eased by taking in consider- 
ation the following: 

a. For small arrays no ambiguities will occur at low 


frequencies. Ambiguities occur when |t|w>r for the largest 


aperture (diameter in circular array). However, for a small 
array when D/A<<l, we get |tlw = ane << m which is unambiguous. 


b. The pairs of sensors have different orientation 
and one can find always (if there are enough sensors) a pair 
of sensors or more which have unambiguous phase difference 
measurement even at higher frequencies than the limit in a. 
Those will be the pairs which are in a line parallel to the 
incoming signal wavefront. 

Using the above two facts the following procedure can 
be outlined to effectively unwrap the phase. 

- Beginning at the low frequencies where there is no 

ambiguity find those pair/s which are perpendicular 


to the bearing of the target/s. 
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- At the higher frequency beginning from those unam- 
biguous pairs, and knowing the geometry, i.e. 
knowing what is the biggest possible change of 
phase between adjacent (in orientation) pairs, 
one can readily detect any 2m jump in phase 
difference and thus effectively unwrap the phase. 
This procedure will be valid only if the target will 
have spectral components under the ambiguity frequency for 
the given array. For example, for a 5 meter diameter array 
this ambiguity frequency is 150Hz and practically we can 
always find some line components below this limit (50 and 


60Hz for instance). 


G. TIME DELAY ESTIMATION 
The time delay estimation process is essential to the 
derivation of the time delay trace function. In some sense 
this process can be more powerful in a detection stage when 
the signals are relatively weak (low SNR per frequency bin) 
because unlike the phase difference trace function it 
"integrates" over all "qualified" frequencies. 
This subject is extensively treated in the current 
" literature and two basically different ways emerged: 
- The Generalized Correlation method which was treated 
in depth by G. C. Carter [6], several works derived and 
extended from it [19] as well as other works on optimal 
and suboptimal versions of correlation schemes [8, 15]. 


- Minimization of some weighted cost function [9, 14]. 
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Both the above methods are closely related to the spectral 
estimation which was reviewed in a previous section (V.E). 

During this work both methods were analyzed, however, it 
seems that the second method is more applicable to the problem 
of small arrays mainly because its time delay resolution is 
not limited by the sampling interval as in the correlation 
method. (In the correlation method the problem can also be 
overcome by elaborated interpolation schemes.) 

As a result the second method will be reviewed here, 
specifically the development done by Chan et. al. [9] which 
uses the weighted least square error criteria. This method 
was adopted and will be followed. 


The model of the signal to be processed is represented by: 


A(t) = C + n; Ct) 


(V-47) 
y(t) = as(t-T) + no(t) 

The signal s(t) is uncorrelated with the corrupting noises 
ni and n, which are stationary random processes. 

Assume that the two inputs x(t) and y(t) are passed 
through a cross spectral estimation process, as described 
before. The phase difference of the cross spectra and the 
magnitude squared coherence function are calculated. These 
two parameters will be used for this processor 


^ 


Ad 


= TW FE n 


E) E 0,1,...M*1 (V-48) 


= D لو‎ TT 
w= MET ٠ 


however the 0 and T frequencies will not be used. 
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7 is an error (noise) with the following statistics: 


"Without signal - E is uniformly distributed from -m 
to m 
T? 
Ete} = 0 Var te} - (V-u9) 


MiB stenal for practical cases of K (EDF > 300) 
we can assume it is approximately gaussian and has a variance 
BT (V-50) 


n p K 


where 


2 
| ۴۶ ا‎ 
ES SD (V-51 ) 


n 2 
l- | Yzy fp? 


The errors for different frequencies are uncorrelated 
l.e. وع) ظ‎ 0 [one] and Gemaier 16 [ م‎ 


The weighted least squares cost function is: 


M 
58 ) 7 (AQ - To ) (V-52) 
n=l 


where the subscripts of En? were dropped. 
M 


۰ ا 
(Ad - Tu Du, (V-53)‏ ¥- ) 2 222 - 0 


n= 1 
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and U Ao نه‎ 


T = (V-54) 


This estimator has two prominent features: 

- It weighs those measurements which have less 
variance more heavily. 

- It implicitly assumes that A ゅ =0 So it provides 
one fixed point on the straight line Ad=ut, a 
feature that intuitively should reduce the variance 
of the estimate. 


The statistics of the estimate are: 


mar} = + + unbiased 


He) - 





which for a reasonable Y is very small for practical values 
of K and M. 
There are some points to be remembered: 
- The factor Y is not given but calculated from 
TT 
xY 


ees — (V-56) 


b el 
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It can be defined by the detection (qualifying) process 
mentioned in the previous section and as such the 
performance of the estimator may be tailored as de- 
sired. 

M 


The summation over the frequency ( £ ) does not 
n-l 


necessarily include all frequencies, only those which are 


qualified 


should be 


by the detection process and for the others ve 


set to zero. 


This section will not be complete without discussing 


the multiple target problem. For a multiple target environ- 


ment the estimator as defined here will not function properly 


because it will not have one distinct line Ad=wt to estimate. 


In order to apply the estimator some preselection of 


the ۵ measurements should be done. 


A 


suggested 


heuristic way to perform the preselection is 

as follows: 

Beginning from the lowest "qualified" frequency 
assume a "target" line (14) which correponds to 
the Ag. at this point and to the zero frequency 
Point. 


Take the next valid ۵۵ and test if 


T1750, A0, «11 *£ (V-57) 


where & is chosen based on some knowledge of the 
variance of AQ. i.e. the threshold on "mE in 


the qualifying process. 
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If it passed the test it is attributed to the same 
"target", Ti. 
If the test is not passed a new "target" line is 


assumed (15) from Af, and the zero frequency. 


Take the next valid Ad. and perform: 
"e n ur ao 0 


1f passed attribute to Ti 


imenot asest for T: 


Ee Es 


1f passed attribute to T, 


if not assume a new target line T 


When this procedure is finished for all qualified 
frequency bins the data points attributed to each 
"target" line may be counted and for those targets 
which pass a certain limit of valid data points 


the estimator in Eq. (V-54) is applied. 


This procedure must be done for all the pairs in the 


system in order to derive the time delay trace function. 


It should be mentioned however that once some targets 


are detected the test will be performed relative to them first 


and if the "new" data does not fit within the limits then new 


"target" lines are assumed. 
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۲ STEM LEVEL PERFORMANCE 

Once the trace functions are derived the bearing estima- 
tors derived in Chapter III may be applied. The resultant 
bearing along with the frequency and power spectral estimate 
for that frequency can be displayed. In this section some 
system level performance curves will be presented. System 
level means that the output beamwidth as defined in Chapter IV 
MENU ted versus the input SNR. The relatims used to reflect 
the phase difference estimation variance back to the input 


SNR. are Equations (V-35) and (V-46) and their assumptions 


set forth: 
3 Y 
SNR; = 27% 
SNR. 
۱۲۱ = eweg 
1 
Eq. (V-35) m و‎ 
2 11د‎ 
g? = 1-|Y| = Er es 
Agij م إا‎ (SNR, ) 
A 1 K 
ZR A 
(SNR.+1) 
1 
2 2 
SNR. +2SNR.+1-SNR. 
9 1 1 1 
5 2 
SNR. KA 
J+2SNR. 
= "m (V-58) 
SNR; K, 
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K, was chosen at a value which will validate the 
gaussian assumption for v(i,j) as defined in Chapter IV, and 
it was set to 250. 

Figures V-4 through V-12 show these performance curves 
for circular and linear arrays along with the same CRLB pre- 
sented in Chapter IV but also related to the input SNR. 

As expected these curves show the same basic behavior 
as the estimator level performance. 

It should be remembered that the relations used, 
expecially Eq. (V-46) are somewhat restrictive in their appli- 
cation and as a result the system performance curves shown 
may not reflect accurately the real situation. Nevertheless 
they can give a good perspective on the expected performance. 
NOTE: The SNR; is per analysis bandwidth and not total SNR 


at the input. 
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VI. CONCLUSION 


1 SUMMARY OF THE PRESENT WORK 

The aim of the work performed in this dissertation was to 
develop a solution to the problem of signal processing, 
particularly spatial processing, of small arrays at low fre- 
quencies. Special emphasis was given to the circular array 
geometry which was motivated by its wide practical use. 

As presently available methods looked unsatisfactory to 
the author, a new concept called here the "trace function" 
was formulated and developed in this work. 

It is based on the fact that the directional information 
of an incoming plane wave is mainly contained in the received 
Signal phases, or to be more precise, in the relative phase 
differences between received signals at the various elements 
of the array. 

Those phase differences are representative of the time 
delay of the plane wave between the elements of the array. 
The notion is that the phase differences/time delays of the 
Signals received from an array of sensors which has a given 
defined geometry describe a certain "trace function" which is 
constant for that array at a given look direction. If one can 
find a geometry for which the trace function's basic shape is 
not dependent on the look direction, but only some of its 


parameters are, then one has a powerful method to correlate 
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(compare) the trace function of the received signal to a stored 
replica. 

The concept of the trace function was formalized and its 
application to some typical array geometries was demonstrated. 
Furthermore, it was shown that for highly symmetric arrays like 
the circular and linear arrays the trace function reduces toa 
particularly simple form. 

This characteristic was used to derive a relatively simple 
and manageable MMSE. estimator for the bearing of the in- 
coming signal. The estimator is applicable to either narrow 
band signal by use of phase difference trace function or to 
the wide band signal using time delay trace function. 

The performance of the estimator was checked by simulation 
and compared to the CRLB as adapted to this application. 

A brief performance summary is: 

~ The estimator exhibits a threshold phenomena called 

"lock-on" with respect to the variance of the phase 
difference measurements. 

- Within the "lock-on" range the estimator is practically 

unbiased. 

- With respect to beamwidth performance the estimator is 

asymmtotically efficient. 

- It exhibits a saturation point in the number of elements 

at constant diameter or in the diameter at constant 
number of element, beyond which only marginal improve- 


ment is achieved in beamwidth. 
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The estimator was applied to both circular and linear 
arrays and its performance was equally good except for 
differences imposed by the geometry. (symmetry) 

Finally, a system configuration applying trace function 
principles was outlined and the major problem areas caused by 
the specific application were identified, reviewed and some 
suggestions to solutions were made. 

As a final conclusion it can be said that this work pre- 
sented a new way of looking at the available information from 
an array of sensors. It is clear that the use of the above 
concept in a generalized replica correlation model is effective 
for highly symmetrical arrays. As such the circular array is 
a natural one to be treated this way. For this array the 
derivations of the bearing estimator is simple and its per- 
formance appears to be very good. 

However, it should be emphasized that the trace function 
concept is applicable to any array geometry and for both 


narrow and wideband processing schemes. 


B. OTHER POSSIBLE APPLICATIONS 

The trace function concept and the estimators derived in 
this work are not limited in their application to the passive 
sonar case. 

They can be applied to related areas like seismic, and 


ocean wave analysis as well. 
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In addition the principle can be applied to any other 
situation in which phase differences and/or time delays 
Within an array of sensor can be measured. The performance 
curves show that the trace function bearing estimator is 
very tolerant to phase difference measurement noise, even 


0 0 ° < e š 
-40” of standard deviation in the measurements can give 


30 
very good results with reasonable number of elements. This 
suggests that even in areas where F.F.T. techniques for cross 
spectral estimation are still beyond the state-of-the-art 
(like microwave frequencies) available hardware exists for 
phase difference or time of arrival measurements which may 


be sufficient to allow application of the trace function 


Eoncept. 


PROPOSED FUTURE WORK 
Although the trace function principle was developed in 
this work there are a vast number of points to be checked. 
In addition several extensions and applications should be 
investigated. 
Some proposed topics are: 
1. Testing the concept of 3 possible ways: 
a. Computer simulation of a whole system 
(similar to Figure V-1) including model of 
the sea. 
b. Laboratory simulation. 
C. Testing on real data at sea. 
2. Investigation of the weighted MMSE , optimi- 
zation of the weights and their possible adaptive 


derivation. 
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3. Derivation of a minimization algorithm for "non- 
analytic" trace function. 
4. Further formalization and development of the 
phase preserving spatial decorrelation. 
5. Investigation of the multiple target discrimina- 
tion for both narrow and wide band targets. 
The above are only a few of the possible topics opened up 


by the introduction of the trace function principle. 
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APPENDIX A 
TRACE FUNCTION DERIVATIONS 
IEMPeulap Array 
This derivation is related to Figure II-4. Using Eq. (II-2) 
the following is the interpretation of the parameters for cir- 


cular arrays: 


a = -sindx + cosdy 

一 9 ee 
بسح‎ == EE + = 

ER [ sinó” ix + cos y iy] 


where x, y are unit vectors in x and y directions. 


Inserting in Equation (II-6): 


AT.. = (-sin®x + cos 0y) * [- (sin 1 = sin“ 


is J) x 


Ol 


27ِ DT 17 
+ (cosT e cosy 3) y ] 


[cos 6 (cost i - cos T s ۳ sine(sin 和 rr 1 - sing! at 


OA 


[cos (£T. i - 0) - cos (<7 j - 65] 


١ 
Ol 


- -2$ sin] (itj) - 9] sin[f (i-j)] (A-1) 
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and for phase difference Eq. (II-5): 


۰ _ AC ۵+ cio ۲ ی وم‎ 7 
و د۵‎ = AT; DER m sin [= Crea 9]sin[x (i-j)] (A-2) 


2. Linear Array 


This derivation is related to Figure II-15. The following 
are the various terns: 
« = sindx + cosdy 


O Î xX 
1 


S| 


1 ۲ 91۳۵ in Eq. (II-6): 


BER. = (sin9x + cos0y):[(i - jOx] 


al 


= $ sine(i-j) (A-3) 


and for phase difference Eq. (II-5): 


s 3بس"‎ - : 5 
Bis و‎ = E w, = I (i-j)sin® (A-4) 
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l. Computer Program for Three Dimensional Plotting of the 


Trace Functions for Circular Arrays. 
The following computer program as well as all the other 
programs in this dissertation were developed and executed on 


a Hewlett-Packard 9845T computer. 
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2. Computer Program for Three Dimensional Plotting of the 
Trace Function for Linear Arrays. 


Following is the listing of the computer program for the 


plotting of the trace function for linear arrays. 
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computer Program for Three Dimensional Plotting of the 


Trace Function for Random Arrays. 
Following is the listing of the computer program for the 
plotting of the trace function for random arrays. The same 


program may be used for conformal arrays. 
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APPENDIX C 


SIMULATION AND THE NOISE MODEL 
NoE MIDDEL 


The simulation of the performance was done by generating 
a noisy trace function (Eq. III-1). 


^ 


AQ; 5 fmp(i,j,0) + v(1,3) (0-1) 


and applying the MMSE estimator. 

The trace function mp(1,3,0) is generated simply by 
calculating its value for all "i" and "3" at a given bearing 
"0". The additive noise v(i,j) is the more problematic term 


and its definitions will be discussed next. 


l. Noise Model 

In Ref. [29] it was shown that for high SNR or conversely 
for low SNR with long averaging time (high K) the error in the 
phase difference estimate may be considered approximately as 
a gaussian distributed noise with zero mean and the variance 


given by Eq. (V-35) 


2 
1-|v..| 
ota : | 1j (C-2) 
ej R 
12 A 


at a specific frequency. 
This assumption of gaussian noise was used throughout the 


simulation and to validate it KA was chosen to be 250. 
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This K corresponds, by using Eq. (V-37), to a T (data 
stationarity time) of approximately 85 seconds when assuming 
one hertz frequency resolution in processing. This figure 
is a reasonable one for real data in the ocean. 

The values of coherence square used were in range - aw 
2929985 which corresponds to input SNR of -15 to +21 dB 
(using Eq. (V-46)). 

The next point to be considered is the covariance between 
the various noise terms. In the phase difference estimation 
process each element is used in 2N-1 pairs so that obviously 
there is a certain correlation between these pairs. The 
following model is assumed for the covariance matrix which 


takes in consideration the various pairs of terms of the 


trace function. Table C-l summarizes the above model. 


DESCRIPTION TERM FORM COVARIANCE 
| Pairs without common elements and و۵۵‎ ws 










Pairs with one common element 






Pairs with two common elements 






Pairs with all common elements 


Table C-1 
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The covariance array is basically a four dimensional array 
which takes in consideration all the above pairs of terms. 
۶ 39 ease of perception and calculation this array was 
described as a two dimensional array which is composed of 
submatrices. 

This covariance matrix is of order mee and corresponds 
to the N* noise terms of the trace function. 

An example of such a covariance matrix for a 4 element 
array (16 term trace function) is shown in Figure C-1. 

The method in Ref. [u] was used to obtain correlated 
gaussian random variables having a specified covariance matrix. 
Specifically, defining the above RES covariance matrix as Q 


' then 


Q = UA U (C-3) 


where U - is a matrix whose columns are the eigenvectors of Q 
A- is a diagonal matrix whose elements on the diagonal 
are the eighenvalues of Q 
U transpose 

The next step will be to generate N independent gaussian 
random variables with zero mean and unity variance. 

This operation can be done by several methods, in the 
present case the polar transformation method was used [4] to 
generate a pair of gaussian independent random variables from 
a pair of uniform random variables which are available as a 


built-in function in the computer. 
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To generate the dependent gaussian variables: 
Z (C-L) 


where V - dependent random variable vector 
4 - independent random variable vector 


As a result the covariance of V is: 
BWN} = v A BZ} AU (C-5) 


2 = I - independent, unity variance 


Pyt) =u Pu q (C-6) 


To get the required noise power V will be multiplied by the 


scalar /Np (noise r.m.s. amplitude): 
No= Mp U Mz 2-7) 


2. Simulation 

Once the additive noise was generated it was added to the 
trace function calculated for the circular or linear array 
Eq. (II-8) and (II-10) respectively. 

Next the corresponding MMSE estimator was applied to 
Eq. (III-23) and (III-24). 

The above operations were run 48 times and the mean and 
variance of the resultant bearing estimates were calculated 
and printed out. In addition a plot of the beamwidth versus 


input SNR was done. 
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The CRLB Eq. (IV-3) and (IV-4) respectively are given 
only with the simulation results for comparison purposes. 

Following is a listing of the simulation program for an 
8 element circular array. In order to change the number of 
elements in addition to change of the parameter N the rele- 
vant arrays in the main program have to be redimensioned. 
(line 240) 

In order to change to a linear array, equations have to 
be changed to the applicable ones in the generation of the 
Signals, the estimators and the CRLB (lines 200, 360, 680, 
NE 16/70, 1680, 1710, 1760, 1770, 1800). 

The transformation matrix for the generation of the 


5١ has to be calculated 


correlated random variable (t = U A 
only once for a given number of elements, then this part of 
the program can be skipped (lines 260 to 290). Care should 
be taken to have a data file (TRAN) defined for the storage 
of the transformation matrix. The subroutine "Symar" (line 


2170) for the eigen-analysis can be found in the Hewlett- 


Packard Numerical Analysis Software Package. 
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This section consists of two parts: 
- References which are explictely mentioned in this work. 
- Additional subject related bibliography which is an 


outcome of an extensive bibliographical search done 
during the course of this work. 
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